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Abstract In local helioseismology, numerical simulations of wave propagation 
are useful to model the interaction of solar waves with perturbations to a back- 
ground solar model. However, the solution to the linearised equations of mo- 
tion include convective modes that can swamp the hclioseismic waves we are 
interested in. In this paper, we construct background solar models that are 
stable against convection, by modifying the vertical pressure gradient of Model S 
(Christensen-Dalsgaard et al, 1996, Science, 272, 1286) relinquishing hydro- 
static equilibrium. However, the stabilisation affects the eigenmodes that we 
wish to remain as close to Model S as possible. In a bid to recover the Model S 
eigenmodes, we choose to make additional corrections to the sound speed of 
Model S before stabilisation. No stabilised model can be perfectly solar-like, 
so we present three stabilised models with slightly different eigenmodes. The 
models arc appropriate to study the / and p\ to p^ modes with spherical har- 
monic degrees in the range from 400 to 900. Background model CSM has a 
modified pressure gradient for stabilisation and has eigenfrequencies within 2% 
of Model S. Model CSM_A has an additional 10% increase in sound speed in the 
top 1 Mm resulting in eigenfrequencies within 2% of Model S and cigenfunctions 
that are, in comparison with CSM, closest to those of Model S. Model CSMJ3 
has a 3% decrease in sound speed in the top 5 Mm resulting in eigenfrequencies 
within 1% of Model S and eigenfunctions that are only marginally adversely 
affected. These models are useful to study the interaction of solar waves with 
embedded three-dimensional heterogeneities, such as convective flows and model 
sunspots. We have also calculated the response of the stabilised models to ex- 
citation by random near-surface sources, using simulations of the propagation 
of linear waves. We find that the simulated power spectra of wave motion are 
in good agreement with an observed SOHO/MDI power spectrum. Overall, our 
convectively stabilised background models provide a good basis for quantitative 
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numerical local helioseismology. The models are available for download from 
http : //www.mps .mpg. de/projects/ seismo/NA4/ . 
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1. Introduction 

Numerical simulations are an important tool to study the effects of surface and 
subsurface solar structures (sunspots, flows, etc.) on solar oscillations. Since 
the wave amplitudes are small compared to the unperturbed background, the 
equations of motion can be linearised about a background solar model containing 
the solar structure being studied. One requirement of linear simulations is that 
the medium through which the waves propagate must be stable against convec- 
tion to prevent unstable modes, which grow exponentially and quickly dominate 
the solution. A commonly used approach is to consider pol ytropic background 
mode ls which are convectively stable by construction (e.g. ICallv and BoedanL 



19931 ). However, the Sun is not a polytrope. 

This work is motivated to satisfy the need to have con vectively stable back- 

groun d models with eigenmodes similar to those of Model S ( Christensen-Dalsgaard et a/l | 



19961 ) . We note that Model S is not a perfect model of the Sun, however it has 
the advantage that it has been extensively tested and used in helioseismology. 

This article is divided into the following sections: Section 2 specifies the 
problem: the geometry, the equations of motion, the wave attenuation model, 
boundary conditions, and the condition for stability. Section 3 outlines the 
strategy for constructing the models and measuring the eigenfrcqucncics and 
eigenfunctions. Sections 4 through to 7 give a detailed description and charac- 
terisation of the eigenmodes of each of the background models that we obtain. 
In Section 8 we implement a mod el of random wave excitation in the S emi- 



spectral Linear MHD (SUM) code (|Cameron. Gizon. and Daiffallahl . 120071) and 



compute the azimuthally averaged power spectra for CSM_A and CSM_B. The 
power spectra are compared to an observed power spectrum from the Michelson 
Doppler Imager onbo ard the Solar and Heliospheric Observatory (SOHO/MDI) 



( Sch errer et al\ . 11995!) . We conclude with a short discussion of the models and 



their foreseen uses. 

2. Specifications of the Problem 
2.1. Geometry 

In this work we arc interested in modelling a relatively small portion of the 
Sun near the solar surface which extends from 25 Mm below the surface to 
2.5 Mm above and 145.77 Mm in each of the horizontal directions. We define 
the height [z] to be negative below the surface and pos itive above, with z = 
given by Model S ( Christensen-Dalsgaard et all . Il996l ). The region is large 



enough that we can study high-degree low-order (n < 4) modes. Relative to the 
entire spherical Sun, however, the size of the region is small. Therefore, in the 
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horizontal direction we can use Cartesian geometry, rather than spherical, so 
that the problem may be solved more efficiently in (horizontal) spectral space. 
We retain the spherical treatment in the radial direction. In this approximation, 
the operators of the problem, where a is any scalar field and A is any vector 
field, defined in Section 12.21 are given explicitly by 



Va = d z az + ik x x + \k y y (I) 
V-A = - 1 d z [(z + R Q ) 2 A z }z + ik x A x x + ik y A y y (2) 

(Z + Hq) 

where the horizontal wave vector is given by k = k x x + k y y. We note here that 
z + Rq is equal to the radial distance from the centre of the Sun. 

2.2. Linearised Wave Equation 

We want to solve for waves propagating through a solar background model in 
the absence of a flow or magnetic field. For adiabatic oscillations the ideal hydro- 
dynamic equations line arised about an arbitrary, inhom ogeneous , background, 



can be written as (e.g., iLvnden-Bell and Ostrikerl Il967f) 



P dU = V(cVV ■ £ + € • Vp) - V ■ (p£)gz (3) 

where £(k, z,t) is the displacement vector, and c, p, p, and g < are the 
background sound speed, pressure, density, and gravitational acceleration re- 
spectively. The operators are specified by Equations (TTJ) and ©. Waves in the 
Sun are attenuated by turbulent convection. We model this by implementing 
an attenuation parameter, as described in Section 12.31 into Equation [3] in the 
following way: 

p(dt + 7) 2 £ = V(c 2 pV ■ £ + € • Vp) - V • {p£)gz. (4) 

We have modelled the attenuation so that it operates both on the displace- 
ment and velocity. This assumes that turbulence in the Sun redistributes the 
displacement perturbations throughout the atmosphere without necessarily in- 
volving the macroscopic (obser yable) velocity. This leads us t o use v = (d t + 7)^ 



as the observable velocity as in lCameron. Gizon. and Duvall| ( [ 2 008) 



In this article the SUM code is used to solve Equation ((4]) ( Cameron. Gizon. and Daiffallahl .[ 

l2007t ) for two types of simulations: to propagate wave-packets and to simulate the 
stochastically excited wave field of the Sun. The simulations use 1098 uniformly 
spaced (0.025 Mm) grid points in the vertical direction and 100 modes in each 
of the horizontal directions. 

2.3. Damping Layers and Wave Attenuation 



We retain the boundary conditions of ICameron. Gizon. and Daiffallahl ([20071 ) 



where the box is periodic in the horizontal direction and the top boundary 
condition is a free surface (the Lagrangian pressure perturbation is zero). In 
addition, at the top and bottom boundaries, "sponge" layers are implemented 
that artificially reduce the energy of the waves to minimise reflection. 
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Figure 1. The top panel shows the damping 7/27T (with T(fc) = 0) as a function of z. The 
top damping layer is much stronger than the bottom. From left to right, the dotted line is the 
top of Model S (zt), the short-dashed line is the surface, the asterisks are the height at which 
the random sources, z», are implemented and the long-dashed line is the effective bottom of 
the box (zj,). The bottom panel is a plot of the attenuation [r(fc)/27r] as a function of URq. 



Waves in the Sun are attenuated by granulation and have a finite lifetime. 
We model the frequency full width at half maximum of the /- mode power using 
T(k) = r»(fc/fc t ) 2 ' 2 , where Tj2w = 100 ^Hz and fc* = 9O2/i? (|Gizon and BirchL | 
2002b. The LHS of E quation Q uses (d t + t) 2 £ « (df + 2jd t )£, whereas 



Gizon and Birchl (|2002l ) use (d t +r)d t £ = (df+Td t )t Therefor e, the attenuation 



coeffic ient used in our equation of motion is half of that used in lGizon and Birch 



(|2002h . The full form of the damping, j(k, z), shown in Figure [TJ is given by 



7 (fc, z) _ T{k) / e^+ - 85 Mm)/[0.2B Mm] ^ for q 535 < z < 3 5 Mm 



2tt 4tt \ e - (2+18 - 54 Mm )/[°- 625 Mm l ^Hz for - 25 < z < -20 Mm. 

The top damping layer introduces a frequency dependence to the eigenmodc 
solutions. High frequency waves have significant energy in the vicinity of the top 
damping layer and are affected more than the low frequency waves that have 
less energy at these heights. Any damping layers will affect the eigen- 
frequencies and lifetimes of the mode, but in this case the lifetimes 
are predominantly dictated by T(k). The parameters for the damp- 
ing layers were guesses which were shown to empirically damp the 
reflected waves sufficiently and not noticeably affect the eigenfre- 
quencies or lifetimes of the modes. The damping layer parameters 
are not optimised and other forms have also been found to work 
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(e.g. lHanasoge. Duvall. and Couvi datl . l2007h . By using the boundary 



value problem (BVP) solver in Appendix [B] we find that the difference in the 
eigenfrequencies between having and not having the sponge layers is less than 
0.5% for the /, p\, and P2 modes and a little higher for the P3 and p^ modes 
(see Appendix ICl Figure I2"3"t0 . If we adjust the range of the top damping layer 
to 0.125 Mm < z < 2.5 Mm we see a maximum 0.5% reduction but only for the 
P4-modes at high frequencies (see Appendix [Cj Figure [23f). 



2.4. Initial Background Model 



We begin with Model S as our background model (starting from any other 
standard solar model would also be possible). Model S extends to 0.5 Mm above 
the surface, however our computational domain extends up to 2.5 Mm so that 
the boundary conditions are sufficiently far from the surface. We extend Model S 
above z t = 0.5 Mm in the following way: 

cq(z) = c s (z t ) for z > z t , (5) 
p (z) = Ps (z t )e- (z - Zt)/[0 - 125 Mml for z > z t , (6) 
PQ (z) = p s (z t )e- (z -^^- 15 Mm l for z > z t , (7) 

where the subscript "S" refers to Model S, the subscript "0" is the extended 
model. The denominators in the exponents are the scale heights of the density 
and pressure respectively at Zt- The only requirement for the extension of the 
background was that it should not increase the wavespeed since we aim to damp 
the waves at these heights to minimise reflection. Thus, the sound speed was 
held constant and the pressure and density smoothly extended. The extension 
is not meant to represent a realistic solar chromosphere and at this height the 
waves will be artificially damped to prevent reflection. 



2.5. Conditions for Convective Stability 



We want to simulate perturbations superimposed on a background model as- 
suming that the evolution is linear. Part of Model S, and therefore the extended 
Model S described above, is super-adiabatically stratified and convectively un- 
stable. This instability is a real property of the Sun resulting in modes that, 
in a linear calculation, grow exponentially in time and will eventually dominate 
the solution. Therefore, we stabilise the background model against convection to 
satisfy the condition d z p > c 2 d z p. We do this by altering the pressure gradient. 
The reason for choosing to modify the pressure gradient is that it affects the 
eigen modes of the model less than cha nges to the sound speed and/or den- 
sity (jCameron. Gizon. and Duvalll . l2008h . We set the pressure gradient in the 



stabilised model as: 



d z p ■■ 



max(cod zy Oo, d 2 po) f° r z < —0.15 Mm, 

max(c§d^/9 — e i, d z p ) f° r — 0.15 Mm < z < 0.1 Mm , 

max(cQd z/ 9 ,d 2 po) for 0.1 < z < 0.325 Mm, 

max(0.99c^d z p , d z p ) for z > 0.325 Mm , 
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where ei = 10~ 5 cgs (at the surface this is ps 0.002 c\ d z po). This formulation 
was arrived at by empirically testing the stability of the simulation with small 
values of t\. An additional constraint was that it should also remain stable with 
an em bedded perturbation (e.g. a sunspot as in Cameron. Gizon. and Duvalll . 



20081 ). This was the smallest value that was found to satisfy these conditions. 



The derivatives, here, are evaluated numerically as 

d z po(zi) = po(zi)la.[po(z i+1 )/po(zi- 1 )}/(zi + i - Zi-i) 
and 

dzPo(zi) = p (z i )ln[p (z i+1 )/p (z l ^ 1 )} /(z i+1 - z 4 _i) 

to achieve a greater numerical accuracy. We have tested that this criterion is 
effective in maintaining stability for simulations for up to ten solar days. 

The stabilisation forfeits hydrostatic equilibrium and introduces gravity modes| 
into the solution. The gravity modes all have low frequencies and can easily be 
excluded from any subsequent analyses. The lack of hydrostatic equilibrium is 
likely to be more consequential. There are different formulations of the oscillation 
equations, those that incorporate the assumption of hydrostatic equilibrium and 
those that do not. We stress that everything presented in this article applies to 
the formulation presented in Equation [4] which was derived from the equations 
of continuity, energy and motion, respectively 

d t p' = -v-(pd t s), 

d tP ' = c 2 (d t p' +&(€-V/>))-ft(€-Vp) and 
p{dt + l) 2 t = -Vp' + p'gz, 

(where the primed quantities are the perturbations), without assuming hydro- 
static equilibrium. Also, the implications for seismic reciprocity ( Dahlen and Trom p.j 
1998) have not been explored and may be important. 



3. Strategy Outline 

Now that we have set out the problem, we outline the strategy involved in 
developing the convectively stable background models presented in this article. 
It is described as follows: 

• Begin with Extended Solar Model S. 

• Convectively stabilise it by changing d z p as described in Section l2~5l This 
results in CSM. 

• Compare the eigenfrequencies and eigenfunctions to those of Model S. 

• We find that the eigenfunctions near the surface, where we are most in- 
terested in modelling, are not well matched and the eigenfrequencies have 
increased. 

We are left with the choice to modify the sound speed and/or the density to 
try to correct the eigenmodes. Since modifying the density has a large effect on 
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the /-mode energy density, we choose to change the sound speed only. Empir- 
ically, we found that increasing the sound speed near the surface improves the 
eigenfunctions: 

• Begin with Extended Solar Model S. 

• Increase the sound speed in the top 1 Mm by 10% (Equation [5]). 

• Convectively stabilise the model. This results in CSM_A. 

• Compare the eigcnfrcqucncies and eigenfunctions to those of Model S. 

• We find that the eigenfunctions are a better match with Model S than CSM 
and the eigenfrequencies are only slightly over-estimated. 

We attempt to correct the eigenfrequencies by introducing a small decrease in 
sound speed in the top ~ 5 Mm, which will reduce the overall speed of the waves 
and thus reduce the eigenfrequencies: 

• Begin with Extended Solar Model S. 

• Take the sound speed profile of CSM_A and introduce an additional de- 
crease in the sound speed of 3% in the top 5 Mm (Equation [5]). 

• Convectively stabilise the model. 

• Compare the eigcnfrcqucncies and eigenfunctions to those of Model S. 

• We find the eigenfrequencies are closer to Model S and the eigenfunctions 
are only moderately further from Model S than CSM_A. This results in 
CSM.B. 

For a smooth transition, a Gaussian function was selected for the sound speed 
changes. The particular parameters were determined by trial- and-error of a 
few guesses to empirically evaluate how they further affected the eigenmodes. 
The comparisons to the eigenmodes of Model S were judged by eye. We cal- 
culated the eigenmodes of the models in two ways. The first used the SUM 
numerical simulations (see Appendix [K§ and the second used a BVP solver (see 
Appendix [Bj . 

As a quantitative measure of the difference between eigenfrequencies of Model S| 
and the featured models, we compute the relative difference of the real part of the 
eigenfrequencies (determined from both SUM and the BVP) to the real part of 
the Model S eigenfreque ncies, uj/lus — 1. These particular Mod el S eigenfrequen- 
cies were calculated as in lBirch. Kosovichev. and Duvalll (2004) using a Cartesian 



geometry and constant gravity. For the modes we are interested in, the geometry 
and radially dependent gravity affect the eigenfrequencies by no more than 0.5% 
(see Appendix [C]) . We measure the difference between the eigenfunctions of 
Model S and the stabilised model in two ways. The first, by calculating the 
relative difference in area under Rc^y/p] near the surface between Model S 
and the respective stabilised model. The second, by calculating the difference 
of the height [z p ] of the uppermost peak of Rc[u zv /p] between Model S and the 
respective stabilised model for each eigenmode (Section [7]). 



4. Convectively Stable Model (CSM) 

Figure [5] shows the relative difference between the stabilised pressure gradient 
of CSM and the pressure gradient of Model S, d z p/d z ps — 1, which is as large 
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Figure 2. The relative difference of the pressure gradient between CSM and Model S as a 
function of height, z. 



as 35% near the surface. We now discuss the effect this change in the pressure 
gradient has on the eigenmodes. 

Figure[3]shows Re[w 2A /p], normalised so that J° 2 * \J i\ v z\ 2 + \v x \ 2 )pdz = 1, 
as a function of z for / and p\ to p4 eigenmodes from Model S and CSM (derived 
using both SUM and the BVP). Recall that the depth of our domain allows us 
to study only up to the p^ mode. The horizontal velocity component of the 
eigenfunctions, v xy /~p, was found to have a similar agreement with Model S. 

We observe that the main effect of the stabilisation on the eigenfunctions 
is to decrease the amplitude of Reft^y'/?] near the surface. Since this is where 
the stabilisation has the greatest effect on the pressure gradient, changes to the 
eigenfunctions in this region are not unexpected. 

Figure S] shows the relative difference of the real part of the eigcnfrequcncics, 
to /cos — 1, for each radial order as a function of wavenumber. The quantitative 
average over 400 < kR & < 900 shows that the increase in the eigenfrequencies 
is less than 2%. The increase in /-mode eigenfrequencies compared to Model S 
can be attributed to the treatment of gravity and geometry of the operators 
(see Appendix [C]). The agreement between Model S and each of the convectively 
stable models will be quantified in Section [7J 

Since it is a necessity to modify Model S, and therefore no subsequent model 
will have exactly the same eigenmodes, we attempt to correct the eigenmodes by 
modifying the sound speed. We found a trade-off between having eigenfunctions 
or eigenfrequencies closer to those of Model S. In model CSM_A (Section [5]) we 
attempt to improve the eigenfunctions and in CSM_B we try to improve the 
eigenfrequencies without affecting the eigenfunctions too much (Section [H]). 



SOLA: validation_vl0.tex; 20 January 2013; 11:53; p. 8 



Seismic Models for Numerical Hclioscismology 



0.15 : 

0.10 : 

0.05 

0.00 

-0.05 

-0.10 

-0.15 
0.15 h 



-0.15 : 
0.15 



1.7mHz 2.1mHz 2.7mHz ~ 

f ::| 




SUM 
BVP 

Model S 



2.7mHz 3.2mHz 4.1 mHz v 




3.7mHz 4.3mHz 5.3mHz 




2.2mHz 2.6mHz 3.4mHz 

Pi 



3.2mHz 3.8mHz 4.7mHz 

P3 



-2 -4 -6 -8 -10-12-14 
z (Mm) 



-2 -4 -6 -8 -10-12-14 
z (Mm) 

Figure 3. The z-dependcncc of the real component of v z ^/p for a number of cigenmodes 
of CSM. The eigenfrequencies for the wavenumbers, fci?Q = 270, 500, 750, are specified by 
colour. The modes have been normalised so that v z is real at 0.2 Mm and have equal integrals. 
The dashed curve shows the eigenmodes from the BVP solution, the dotted curve shows the 
eigenmodes from the SUM simulations and the solid curve shows the Model S eigenmodes. 
Each panel corresponds to a different radial order [/, p\ to p^\. 
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Figure 4. The relative difference between the real part of the CSM eigenfrequencies, uj, and 
the Model S eigenfrequencies, wg, as a function of fci?Q. The solid curves are for the simulated 
SUM eigenfrequencies and the dashed curves are for the BVP solutions. The average relative 
difference of each radial order in this range is within 2%, with at most 0.5% due to the different 
treatment of gravity and geometry (see Appendix [Cjl . 



5. Convectively Stable Model A (CSM_A) 

We follow the procedure set out in Section [3] We found that an increase in sound 
speed improved the match between the eigcnfunctions of CSM and Model S near 
the surface. We chose 



ca(z) = Co (2) 



1 + 0.1exp - 



1.0 Mm 



(8) 



where the subscript "A" indicates CSM_A. Starting from Model S with ca 
specifying the sound speed, we then rederived the pressure gradient required for 
stability as set out in Section l2~5l Figure [5] shows the relative difference between 
CSM_A and Model S of the sound speed squared and the pressure gradient as a 
function of height. This change in sound speed was found to raise the height of the 
uppermost peak of Re^^/p]. Figure [5] shows Refi^y/jo] for various eigenmodcs 
from CSM_A for each radial order, / and p\ to p±. Particularly, the /-mode 
eigenfunctions are close to Model S. The pi and P2 modes are also a better 
match, especially near the surface. 

The real parts of the eigenfrequencies, shown in Figure are not significantly 
affected: the average (over 400 < kR & < 900) relative difference for each radial 
order is still less than 2% of Model S values. We have constructed a convectively 
stable model, CSM_A, with eigenfunctions closer to Model S than CSM and 
reasonably similar eigenfrequencies. 
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Figure 5. The relative difference between CSM_A and Model S of (a) the sound speed squared 
and (b) the pressure gradient, as a function of z. 



6. Convectively Stable Model B (CSM_B) 

Starting from Model S and ca, we constructed a model with eigenfrequencies 
closer to Model S than CSM_A and reasonable eigenfunctions (as described in 
Section [3]). The eigenfrequencies are related to the phase speed of the wave (w/fc) 
and so we slowed the waves down by adding a broad reduction in sound speed 
of CSM.A. We chose 



c B (z) = c A (z) 



1 -0.03exp 



5.0 Mm 



(9) 



where subscript "B" indicates CSM_B. Figure [5] shows the relative difference 
between CSM_B and Model S (a) sound speed squared and (b) pressure gradient 
as a function of height. 

The eigenfunctions arc slightly adversely affected as can be seen by comparing 
Figure [5] with Figure[51 however they are still more solar-like than those of CSM 
(Figure [3]) . The real parts of the eigenfrequencies (Figure [TU]) reduce to within 
1% of Model S. We have not found a model which resulted in more similar 
eigenfrequencies without grossly changing the eigenfunctions. With this sound 
speed profile, we have arrived at a convectively stable model, CSM_B, with 
eigenfrequencies closer to those of Model S than CSM or CSM_A. 
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Figure 6. The z-dependence of the real component of v Z sJ~p for a number of eigenmodes of 
CSM_A. The eigenfrequencies for the wavenumbers, fci?Q = 270, 500, 750, are specified by 
colour. The modes have been normalised so that v z is real at z = 0.2 Mm and have equal 
integrals. The dashed curve shows the eigenmodes from the BVP solution, the dotted curve 
shows the eigenmodes from the SUM simulations and the solid curve shows the Model S 
eigenmodes. Each panel corresponds to a different radial order [/, pi to pn\. 
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Figure 7. Left: the relative difference between the real part of the CSM_A eigenfrequencies, 
oja and the real part of the Model S eigenfrequencies, ljq as a function of IcRq. The solid curves 
arc the differences using the eigenfrequencies calculated from SUM and the dashed curves from 
the BVP. Right: the relative difference between the real part of the CSM_A eigenfrequencies, 
WA, and the real part of the CSM eigenfrequencies, lj, calculated by SUM as a function of fci?Q 
brought about by the increase in sound speed. 
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Figure 8. The relative difference between CSM_B and Model S of (a) the sound speed squared 
and (b) the pressure gradient, as a function of z. 



7. Comparison of Eigenfunctions 

Quantitatively, we compare the eigenfunctions by finding the relative difference 
of the area under Re^zy 7 ^] between Model S and each convectively stable back- 
ground in the near-surface layers, —1.0 Mm < z < 0.5 Mm. The difference is 
defined by 

p_ ffT(g^M -Rc[v zS ^]) 2 dz 

U f 0.5 Mm,,, r n2 , ■ V iU J 

J-i Mm (Re^sVPsD 2 dz 

This integration range was chosen because this is where the stabilisation has 
greatest effect. For each radial order we take the mean of D over 400 < kR & < 
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Figure 10. Left: the relative difference between the real part of the CSM_B eigenfrequcncics, 
u>b, and the real part of the Model S eigcnfrcquencies, tog as a function of URq. The solid 
curves are differences in the eigcnfrcquencies calculated using SUM simulations and the dashed 
curves are from the BVP. Right: the relative frequency difference of the real part of the 
eigenfrequcncics, calculated by SLiM as a function of wavenumber between CSM_B and CSM_A 
brought about by the reduction in sound speed. 



900, giving a quantitative measure of the differences between the eigenfunctions 
of Model S and the stabilised model [(D)]. Figure [TT1 shows that for the /, p\ 
and P2 modes CSM_A (triangle) has eigenfunctions closest to that of Model S, 
while CSM (asterisk) has those farthest from Model S. 

In addition, we measure the height of the uppermost peak, z p . of Re[w zv /p]. 
Figure [T^] shows z p for each radial order and each model as indicated. From 
this we see that stabilising the background causes z p to drop in height (i.e. the 
difference between the solid and the long-dash curves). Increasing the sound 
speed in a narrow region close to the surface (CSM_A) pushes the peak back 
towards the surface (dotted curves). The broad decrease in sound speed added 
in CSM_B does not change the location of the peak too much (short-dash curves) . 
The sudden transition to very high upper turning points at high wavenumber 
for Model S (particularly for the p\ and P2 modes) is due to the protuberance in 
the Model S eigenfunctions very close to the surface (for example, the / and pi- 
modes in Figure[3]) which is absent in the stable models. The protuberance is due 
to rapid changes in the density scale height close to the surface that disappears 
after the stabilisation (reduction of d z p). 

We now have three convectively stable solar models each having similar, but 
slightly different, eigenfrequencies and eigenfunctions to Model S. Having models 
focused on achieving slight variations of the same goal (more similar eigenfunc- 
tions or eigenfrequencies) gives us the possibility of testing the sensitivity of 
helioseismic analysis techniques to the background properties. 
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Figure 11. The relative difference in area under Re[ti z ^/p] (see Figures I3I6I9H averaged over 
400 < kRo < 900, (D), between each background model - CSM (asterisk), CSM.A (diamond) 
and CSM_B (triangle) - and Model S as a function of radial order. 



8. Modelling the Random Wave Field 

8.1. Random wave excitation model 

We model the random wave excitation by imposing a vertical force, f z , to the 
right-hand-side of Equation The force is specified by 



l /d 2 



(11) 



where k; is a horizontal wavevector, ojj is an angular frequency, d = 0.075 Mm 
is the width of the source, and the acceleration Gij is a realisation of a complex 
Gaussian random variable with zero mean and variance E[\Gij | 2 ] = exp [— (w J ) 2 /2(T 2 ]| 
where <j/2tt = 2.12 mHz (jGizon and Birchl . 120041 ) . The height of the sources is 
at = —0.75 Mm, which is close to the high ly superadiabatic layer where s olar 
waves are expected to be strongly excited ( Nigam and Kosovichev . 1999). In 
reality, the sources in the Sun will also have a wavenumber dependence which 
we have not included. In practice, the sources are generated before the simulation 
commences and saved with a 30 second cadence. The forcing is applied at each 
time step (in cases herein this is approximately 0.13 solar seconds), with the 
value of the applied forcing changing every 30 solar seconds. We remark that we 
first tried to use a Lor entzian for the frequency dependence ([Title et all Il989i : 
Gizon and Birchl . 120021) , corresponding to sources which decay exponentially in 
time. We found that the resulting power was too strong at high frequencies 
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Figure 12. The height of the uppermost peak of Kc[v zy /p\ (see Figures 1.316191 for each radial 
order as a function of kR^. The stabilisation reduces the height of the peak from Model S 
(long-dash) to CSM (solid). The consequence of adjusting the sound speed is shown in CSM_A 
(dot) and CSM.B (short-dash). 



compared with observations, and that the Gaussian distribution produced a 
better agreement. 



8.2. Azimuthally Averaged Power Spectra 



In this section we used SUM to investigate the response of CSM_A and CSM_B 
to the random wave excitation model as described in Section [8TTI A total of 16 
hours was simulated, however the first eight hours, during which the wave field 
is reaching a steady state, are discarded. To mimic SOHO/MDI observations, we 
save vertical- velocity data at a h eight o f 0.2 M m above the surface (the height at 
which SOHO/MDI observes, see lBruld (|l993l) 1 and account for the modulation 
transfer function of the instrument by m ultiplying the simulated power spectra 
by th e modulation transfer function of Rabello-Soares. Korzennik. and Schoul 
<|200lh . 
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To make a comparison with an observed power spectra, we took eight hours of 
Postel projected (centred at a longitude of 170° and latitude of —8.3°) full-disk 
Doppler observations with a 60 second cadence from SOHO/MDI on 21 January 
2002. The observations consist primarily of quiet Sun covering a surface area 
identical to the simulations. 

We consider the azimuthally averaged (with bin size Ak = 27r/[145.77 Mm]) 
power spectra of the observations, P(k x ,k y ,uj) = \v\ os (k x , k y , w)| 2 , and of the 
simulations with CSM_A and CSM_B, P(k. x ,k y ,uj) = \v z (k. x , k y , lu)\ 2 are shown 
in Figures [T2HH1 and IT51 respectively. The dashed curves are the cigcnfrcquencies 
calculated from Model S for comparison. The straight solid line is where uj/k is 
equal to c(z(,)/(l + zi,/Rq) and zj, = —22.6 Mm; as stated previously, modelling 
a higher co/k would require a deeper box. There is some power evident in the 
low frequencies which are most likely (/-modes introduced by stabilising the 
background. These are the artificial product of having a stable model. Thus, this 
region cannot be compared to solar observations. The remaining "comparable 
domain": b(k) < ui/2ir < k c(zb) j (1 + Zf, / Rq) where b(k) is the lower curve shown 
in these figures, should contain modes which are comparable to those on the 
Sun. The azimuthally averaged power spectra are normalised to the mean power 
within a region defined by {kR - 600) 2 /200 2 + (w/2tt- 3 mHz) 2 /(l mHz) 2 < 1. 
By inspection, the power spectra of CSM_A (Figure PH)) and CSM_B fFigurc[T5|) 
look qualitatively similar to the observed spectrum (Figure \T3§. We now take a 
closer look at the properties. 

8.3. Amplitudes of the Power Spectra 

Figure [If)] shows vertical cuts through the power spectra in Figures Q21 HH and 
1151 as a function of frequency The Model S eigenfrequencies (vertical lines) are 
larger than those of the observations (solid curve), while CSM_A (dash curve) 
and CSM_B (dot curve) eigenfrequencies are larger than those of Model S. It 
also shows that the maximum power and linewidths of the ridges agree with 
observations best at low frequency. 

Figure IT71 shows the total power in the comparable range for the observations 
(solid curve), CSM_A (dash curve) and CSM_B (dot curve) as a function of (a) 
frequency and (b) kR@. The maximum power in the simulations occurs at a 
larger wavenumber than in the observational power. Correcting this could be 
done by fine tuning the wave excitation model, and may be done in the future, 
however the results presented here are sufficiently close for a large number of 
studies. 



8.4. Fitting the Power Spectra 



We analyse the properties of the azimuthally averaged po wer spectra in Fig 
ures 1131 1141 and 1151 by fitting asymmetric Lorentzians (e.g. 
GizorIT2006l) . 



Duvall et al.. 1993 



L(u) = J>« 



n=0 



(1 + B n X n ) 2 + B 2 n 

1 + X2 



N 



(12) 
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Figure 13. The azimuthally averaged power spectrum of eight hours of quiet-Sun SOHO/MDI 
Dopplcr observations. The eigenfrcqucncics of Model S are the dashed curves. The straight 
solid line and the bottom solid curve form the boundaries of the comparable domain of the 
simulations. Stronger power is black and weaker power is white. 



where X n = — - and B n = 



r„/2 u> n — lu v 

to cuts at fixed wavenumber as a function of frequency. In Equation (|12[). the 
maximum power of the rfi 1 ridge is given by P n and is located at a frequency lu„ , 
the valley is at u v , the noise is N, and the full- width-at- half- maximum (FWHM) 
of the asymmetric Lorcntzian is T n . The fitting is done using a Levcnberg- 
Marquardt algorithm for least squares curve fitting using the IDL mpf it package. 
The frequency range of the fit is from w 0.6 of the /-mode Model S eigen- 
frequency to w 1.1 of the p4-mode Model S eigenfrequ ency. We define the 
asymmetry parameter as Xn = B n u n /{T n /2) ( Gizon . 20061) . 
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Figure 14. The azimuthally averaged power spectrum of eight hours of simulated random 
wave excitation in CSM_A, accounting for the SOHO/MDI modulation transfer function and 
presented on the same log-power scale as Figure [T3l The eigenfrequcncies of Model S are the 
dashed curves. The straight solid line is where ui/k is equal to c^(zi,)/(l + z^/Rq), and the 
bottom solid curve [b(k)] form the boundaries of the comparable domain. 



Figure [18] shows the maximum power of each n from fitting Equation (|12[) 
to the power spectrum of the observations (top), CSM_A (middle) and CSM_B 
(bottom). The simulated power spectra have stronger power at high frequency 
than the observations. In addition, the maximum power of n = 1 occurs at a 
lower frequency in the simulations than in the observations. 

Figure \W\ shows the FWHM of the Lorentzian fit for each mode in the power 
spectrum of the observations (top), CSM_A (middle) and CSM_B (bottom). 
The FWHM of the ri dges in the observations is consistent with Figure 2 in 
Antia and Basu (Il999h . keeping in mind that these are coarse measurements. 
The simulation ridges have larger FWHMs than the observations for / and pi 
modes. 
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Figure 15. The azimuthally averaged power spectrum of eight hours of simulated random 
wave excitation in CSM_B, accounting for the SOHO/MDI modulation transfer function and 
presented on the same log-power scale as Figures 1131 and 1141 The eigenfrequencies of Model S 
are the dashed curves. The straight solid line is where ui/k is equal to cg(zi,)/(l + z^/Rq), 
and the bottom solid curve [&(&)] form the boundaries of the comparable domain. 



Figure [501 shows the relative difference of the central ridge frequencies to 
Model S for the observations (top), CSM_A (middle) and CSM_B (bottom). The 
results from the Lorentzian fitting are within 1% of the BVP solutions as shown 
in Figure |2"T1 

Figure [22] shows the \n asymmetries of the observations (top), CSM_A (mid- 
dle) and CSM_B (bottom). We achieve the correct sign and comparable magni- 
tude of the asymmetry for all the modes. The /-mode has negative asymmetries, 
and the value of the asymmetries increases with increasing mode number which 
is in agreement with Gizon ( 20061) . 

We have demonstrated the response of the numerical simulations of wave 
excitation in the Sun using two of the convectively stable background models, 
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Figure 16. Cuts (smoothed over 0.035 mHz for the purpose of this plot) through the az- 
imuthally averaged power spectra in arbitrary units at the indicated wavenumbers for CSM_A 
(dash), CSM_B (dot) and observations (solid) as a function of frequency. The vertical black 
lines are the eigenfrcquencics of Model S. 



CSM_A and CSM_B. The eigenmodes of the background models and the param- 
eters of the sources of acoustic wave oscillations are sufficient to be used as a 
foundation for quantitative solar-like simulations. 

In addition, we have successfully implemented the stable background models 
into the framework of another code which also computes linear simulations of 

helioseismic wave propagation, the Seismic Propagation through Active Regions 

and C onvection (SPARC) code ( Hanasoee et al. . 2006 : Hanasoee. Duvall. and Couvidati .[ 



2007). 



9. Discussion 



We have created three convectively stable solar models which, to slightly differing 
extents, have similar eigenmodes to those of Model S. We have also computed 
helioseismic simulations using a model for the random excitation of waves, which 
together with the stable solar models, reproduce the SOHO/MDI observed mode 
frequencies and asymmetries well for each of the / and p\ to p± ridges. The 
linewidths of the ridges and the power distribution are reasonably similar to 
those of the Sun. 

Although stabilising the background model is an importa nt step in numerical 

studies of wave propagation (and has been done before, e.g. bv lParchevskv and Kosovichevl [ 
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Figure 17. The total of the azimuthally averaged power in the comparable range as a function 
of (a) frequency (averaged over wavenumber in the comparable range) and (b) IcRq 
(averaged over all frequency in the comparable range) for CSM_A (dash), CSM_B 
(dot) and observations (solid), in arbitrary units. 



2007; Cameron. Gizon. and Duval! , 20081: Shelvag. Fedun. and Erdelvi , 2008 : Schunker. Cameron, and Gizon 



2010I ). its effects on the eigenfunctions and eigenfrequencies has received little 
attention. An optimal way to produce a convectively stable background model for 
numerical simulations has not been formulated, but nevertheless the models pre- 
sented here should be useful for a range of studies. In particular, we envisage that 
these models will be used to study the propagation of solar waves through three- 
dimensional h eterogeneities, s uch as convective flows, granulation and model 
sunspots (e.g. Cameron et aZ.I . [201 1 : Dombroski. Birch, and Braun . 2011 ). Hav- 
ing three models with slightly different properties will enable us to quantitatively 
test the sensitivity of the results to the details of the models. The models and 
extra information from the analysis in this paper arc available for download from 
the HELAS local helioseism ology website ( http://www.mps.mpg.de/projects/seismo/NA4/il 
Schunker and Gizonl . l2008h . 
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Figure 18. The maximum power, P n for n = 0, 1, 2, 3 ridges calculated by fitting Equa- 
tion d 1 211 to the azimuthally averaged power spectra as a function of frequency. The top panel 
shows results from the observations, the middle panel from CSM_A and the bottom panel from 
CSM_B. Each ridge is presented by a different symbol as indicated in the legend. 
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Figure 19. The FWHM, r„, for n = 0, 1, 2, 3 ridges calculated by fitting Equation ([T2J to the 
azimuthally averaged power spectra as a function of frequency. The top panel shows results 
from the observations, the middle panel from CSM_A and the bottom panel from CSM_B. The 
symbol legend is the same as in Figure [TBI 
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Figure 20. The relative difference of the central ridge frequencies [ui n ] for n = 0, 1, 2, 3 ridges 
calculated by fitting Equation (1 1 2 1> to the azimuthally averaged power spectra, to those of 
Model S as a function of fci?Q. The symbol legend is the same as in Figure [181 The top panel 
shows results from the observations, the middle panel from CSM_A, and the bottom panel 
from CSM.B. 
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Figure 21. The relative difference between the eigcnfrcqucncics of the BVP solutions, /bvp> 
and the frequency of the maximum ridge power as identified from fitting the power spectrum, 
/fit i for CSM_A. The symbol legend is the same as in Figure fTSl 
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Figure 22. The \ n asymmetries for n = 0, 1, 2, 3 ridges calculated by fitting Equation (1 1 2 1 > 
to the azimuthally averaged power spectra as a function of fci?0. The symbol legend is the 
same as in Figure 1181 The top panel shows results from the observations, the middle panel 
from CSM.A and the bottom panel from CSM.B. 
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Appendix 

A. Calculating Eigenmodes Using Simulations 

The following procedure calculates the eigenmodes of the system using SUM and 
is designed to be applied iteratively. We began by simulating the response of the 
system to a wave packet constructed from Model S eigenmodes of one radial 
order (as in Cameron, Gizon, and Duvall, 2008). The outputs, v x (k,z, t) and 
v z {k, z, t), of a five hour long simulation were saved with a one minute cadence. 
We then took the Fourier transform of the velocity field in time, v z (k, z,w). 
From this we determined the eigenfrequencies from a linear fit in time to the 
phase <j){k, t) = Aig[v z (k, z — 200 km, t)] with the 2ir wrap-around removed. The 
function we used to fit the phase, 4>(k, t), is given by (p(k, t) = Re[L0i(k)]t+(j) o g (k). 
The h e ight o f 200 km corresponds to the observation height of SOHO/MDI 
(|BrulsL 1 19931 ). We determined the radial component v z {(k, z,uj) by applying a 



broad ridge filter to isolate the appropriate radial order, centred on the improved 
(subscript i) estimate of the real part of the eigenfrequencies Re[wi(fc)]. The same 
filter was applied to the horizontal velocity to get v x f(k, z, lo). The velocities are 
then Fourier transformed from frequency space back to time. 
The improved eigenmodes are then given by 

J^ 0amm v xf (k,z,t)exp[-i(uj i (k)t + (j) oa )}dt 

Vxi(K,Z) — 300min 



v zi (k,z) = 



Jo mm VyA z = km, t) cxp[—i(oji(k)t + <j) g)]dt 

Jo mm v z t(k,z,t) exp[-i(u)j(k)t + (j) ff)}dt 
j mm Vzf (fc^ z = 200 km, t) exp[— i(u>i{k)t + (j) e)]dt 



Note that (f> s{k) and the denominator arc defined from the vertical velocity 
component, v z i(k, z) = 1 at z = 200 km. From these eigenmodes we constructed 
a new wave packet initial condition and the simulation was re-computed with 
this wave packet. In practice we found that a single pass is sufficient and the 
improved eigenmodes from the first simulation were used to compare to Model S. 



B. Determining the Eigenmodes of the Boundary Value Problem 



The perturbations of a particular eigenmode with radial order n of the Equa- 
tion ((I]) have the form 



v(k,z,t) = [v z (k,z)z+v x (k,z)k]e- i ^-' k -^ 
j>'CM,i) = P '(k,z) e - l ^ t - k - x > 



(13) 
(14) 



with v„ = (d t +7)£„- 

After some manipulation (using the continuity equation, equation of motion 
and energy equation) , our system of equations becomes 



pf3v z = - 



dp' 

d7 



pv z v z dp p dv z k 2 p' 



(15) 
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2pv z pv z d7 p dv z p'k 2 
~rf ~ Jpdz + /3~cL7 + 7^ 



v^dp 
j3 dz' 



(16) 



where j3 = 7 — ko. 

Following the method ol lBirch. Kosovichev. and Duvalll (2004), we substitute 



ip 

2/1 = —p= 

2/2 = Vzy/pc 



into Equations (fT5|) and ([T6")) to get 

z 2 k 2 



yivpc 



1 



1 



P 2 



111 



pc \ r/3 



2c 2 p c 2 p d-y 1 dp 
~ ~W~dz~ ~ /3dz" 



c 2 p d 



2/2 



(17) 



and 



gk 2 yisfpt 

/3 2 i 



2/2 



/pc 

5P_d_ 

' dz 



r/3 



in 



(3dz 
d_ 1 y 
dz 



P d7 



= 0. 



(18) 



Then multiplying Equation (jTTJ) by ^fpc and Equation (|18[) by i/y/pc and rear- 
ranging, we get 



dj/2 / 1 d7 1 1 1 dp 
~d7 ~ V2 \J3~dz ~ 2H C ~ 2H~ P ~ ~p~^dz 



and 



dj/i 
dz 



-2/1 



1 

2i^ 



2H„ 



12/2 



r/3c 



g dp 

/3c 3 p dz 



/3 c/fc 2 

- + T 



9 dp 

/3pc dz 







(19) 



(20) 



where 1/H C = —d z c/c and l/ iif = —d z p/p. Equation ([TO]) and (12"U1) reduce to 
Equations (A10) and (All) in lBirch. Kosovichev. and Duvalll (|2004l) in the case 
where the attenuation is not dependent on z, the background is in hydrostatic 
equilibrium and the geometry is Cartesian. 

The top boundary condition is a free surface such that the Lagrangian pressure 
perturbation [Sp] is zero. This means that p' = — £ ■ Vp. The bottom boundary 
is specified by v z — and p' = 1. The boundary conditions translated to 2/1 and 
2/2 are that pcy\ + i?/2/3d 2 p = at the top and 2/2 = and y\ = 1 at the bottom. 

We solve this boundary value problem using the Matlab program bvp4c. In or- 
der to be consistent with the eigenfunction solutions from the SUM simulations, 
we do a similar normalisation of the cigenfunctions so that V z {k, z = 200 km) = 
1. 
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C. Solutions to the BVP for Different Background Models 

We use the BVP solver outlined in Appendix [B] to explore the effects on the 
eigenfrequencies by changing different parameters of the problem with CSM_A. 
To test the robustness of the BVP solver we added 1% noise to the eigenfrcquency 
guess that results in a relative difference of less than 10 -5 as shown in Figure 1251 
(a). In Figure [23] (b) we do not apply any wave attenuation, i.e. 7 = 0. The 
eigenfrequencies decrease in value compared to the CSM_A eigenfrequencies, 
more so for the higher order modes. In Figurc[23j(c) we have set a constant grav- 
itational acceleration of g = — 273.98m/s 2 . This mostly affects the /-mode, but 
the eigenfrequencies are also decreased for the p-modes. Removing the sponge 
layers, so that 7 = T(k), give results, Figure 1251 (d), that are similar to (b). 
Using the full Cartesian operators, as opposed to the spherical derivative in the 
radial direction as in Equation[2j affects the eigenfrequencies at low-wavenumber 
the greatest, as shown in Figure [251 (c). In Figure [251 (f) we have lowered the 
top damping layer to have ~/{k,z)/2ir = r(fc)/47r + e [( z+1 - 28 Mm )/0-25 Mm]^ for 
0.125 < z < 2.5 Mm (retaining the bottom damping layer), which decreases the 
eigenfrequencies. These frequency shifts are small compared to the frequency 
shifts caused by the convectively stabilising the models. 
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Figure 23. The relative difference between the CSM_A eigenfrequencics with various modified 
quantities, uj q , and the BVP eigenfrequencics of the original CSM_A, u. The panels show the 
relative difference for with (a) 1% noise added to the eigenfrequency guess, (b) no damping 
layers or attenuation, (c) constant gravity, (d) no damping layers, but retaining the attenuation, 
(e) Cartesian geometry and (f) top sponge extended lower in height. The frequency shifts are 
much smaller than those introduced by the convectively stable models. 
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